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Description 

FIELD OF THE INVENTION 

5 [0001] This invention relates to Volume Graphics and more specifically to a discrete object representation that allows 
accurate reconstruction of surface normals for Volume Rendering, haptic exploration of volumetric models using force 
feedback, and physically-based calculations of interactions between graphical objects. 

BACKGROUND OF THE INVENTION 

10 

[0002] In Volume Graphics, objects are represented as arrays of sampled data elements. Volumetric object represen- 
tations have some advantages over surface-based graphical object representations because they can represent object 
interiors. This is important both for visualizing the internal structure of complex objects and for physically-based mode- 
ling of actions such as object deformation or the cutting and tearing of volumetric objects. Using volumetric object rep- 

is resentations to model physically-based object interactions is discussed, for example, by S. Gibson in "Beyond volume 
rendering: visualization, haptic exploration, and physical modeling of voxel-based objects", in Proc. Eurographics work- 
shop in Visualization in Scientific Computing, pp. 10-23, Chia, Italy, 1995. Using a volumetric object representation to 
model object deformation is discussed by S. Gibson in "3D ChainMail: a fast algorithm for deforming volumetric 
objects", in Proc. Symposium on Interactive 3D Graphics, Providence, Rl. pp. 149-154, 1997. 

20 [0003] Volume Rendering is technique used for displaying volumetric objects on a 2D computer monitor. There are 
several approaches to Volume Rendering, including ray casting where rays from each pixel in the image plane are cast 
into the object accumulating color as they pass through the volumetric object; and splatting where each element in the 
object is projected onto the image plane, contributing color to the rendered image. For Volume Rendering, information 
about the color, transparency, surface normal, and material reflectance properties for each volume element are needed. 

25 These are often derived from two values stored for each element: a measured sample intensity or density value, and a 
material classification. For example, the density value might be used to calculate the element's opacity and surface nor- 
mal, while the classification might be used to assign colors and reflectance properties to object elements. 
[0004] In computer graphics, the quality of a rendered image is greatly improved by adding an illumination model to 
the rendering process. An illumination model allows the calculation of lighting effects such as spot lighting from specific 

30 light sources, reflections from object surfaces, and shadows. In computer graphics, a typical illumination model is the 
Phong illumination model, which expresses the light reflected from an object surface as follows: 

lr-KJJ d- *K - IJ^*K i l L l cN (PVJ n i 

35 where l n l fl , l^, and l L are the reflected, ambient, object, and light source colors, K a , Kd, and Ks are the ambient, diffuse, 
and specular reflection coefficients, N is the surface normal of the object, L is the direction of the light source, V is the 
viewing direction, and R is a function of N and L 

[0005] In the Phong illumination model, an accurate representation of the object surface normal is essential. Inaccu- 
rate surface normals can have a large effect on the perceived shape and texture of a graphical object. Since, in most 
40 graphical representations, the surface normal is known accurately at only a relatively small number of positions, esti- 
mating the surface normal at arbitrary object positions is very important for high quality rendering. In Volume Graphics, 
estimating surface normals is especially challenging because accurate surface normals are not usually stored in the 
data. Instead, surface normals are usually estimated from data intensity or density values. 

[0006] It is important to note that the calculation of accurate surface normals in Volume Graphics is not only needed 
45 for Volume Rendering. An accurate surface normal estimation is essential for haptic data exploration where a force 
feedback device is used by the user to factually explore the surface of a graphical object. Accurate surface normal esti- 
mation is essential in this application because the sense of touch is even more sensitive than the eye to surface irreg- 
ularities. An accurate surface normal is also important in modeling physically-based interactions between objects 
because the calculation of impact and frictional forces between objects depends on these normals. 
so [0007] Many methods have been used to add lighting effects, or shading, to Volume Rendered images. These meth- 
ods depend on the data representation. Different methods are used when the data is represented as a binary field 
where data samples are either empty or full of a given material, or as a sampled continuous field, such as a sampled 
intensity or material density field. Many methods for shading volumetric objects are discussed, for example, in "Volume 
Visualization", edited by Arie Kaufman, and published by the IEEE Computer Society Press, Los Alamos, CA, in 1991 . 
55 The more sophisticated of these methods require surface normal estimates at sampled data points. 

[0008] In order to estimate surface normals when the object is represented as sampled continuous field, Volume Ren- 
dering often takes advantage of the fact that the image intensity or material density values change fastest at object sur- 
faces. The directional change in image intensity can be represented by the 3D image gradient vector, which, at any 
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image position, points in the direction of fastest change in intensity values and which has a magnitude proportional to 
surface strength. Hence, this gradient vector is often used in Volume Rendering to estimate surface normals for illumi- 
nation. 

[0009] There are many ways to estimate image gradient vectors from sampled data. One straight-forward method is 
5 called the central difference method. Here the gradient vector, (g x , g y g 2 .) is estimated from the 3D image intensity I (x, 
y. z) by: 

g x = l(x+i,y,z)-l(x-i,y,z), 
w g y = l(x,y+1,z)-l(x,y-1,z), 

g 2 = l(x,y,z+1)-l(x,y,z-1). 

[001 0] Other image gradient operators include Sobel operators, cosine filters, and more complex filters such as those 
is discussed by S. Marschener and R. Lobb in "An evaluation of reconstruction filters for Volume Rendering", Proc. IEEE 
Visualization *94, pp. 100-107, Washington, DC, 1994, by T. Moller, R. Machiraju. K. Mueller, and R. Yagel in "Evaluation 
and design of fitters using a taylor series expansion", IEEE Transactions on Visualization and Computer Graphics, Vol. 
3, No. 2, pp. 184-199, 1997, and by M. Bentum, B. Lichtenbelt, and T. Malzbender in "Frequency analysis of gradient 
estimators in volume rendering", IEEE Transactions on Visualization and Computer Graphics, Vol. 2, No. 3. pp. 242- 
20 254, 1996. 

[001 1 ] S. Wang and A. Kaufman calculate smooth surface normals from volumetric data in "Volume-sampled 3D mod- 
eling", in IEEE Computer Graphics and Applications, pp. 26-32, Sept. 1 994. They start with binary segmented data and 
the apply a Gaussian filter to the object intensity or density values. They then estimate surface normals as the gradient 
of this smoothed image. 

25 

SUMMARY OF THE INVENTION 

[0012] In contradistinction to the above methods, in the subject invention, a volumetric distance map in which dis- 
tances to the nearest object surface are stored for each volume element or voxel, is used to improve the surface normal 
30 estimation for volumetric objects. After having created a discrete volumetric distance map, one can calculate a gradient 
vector of the distance field for any point in the volume. Because the gradient vector derived from this distance map 
points in the direction of the surface, it provides the surface normal direction. The distance itself can be used to derive 
the surface normal magnitude. 

[0013] While others have used volumetric distance maps in Volume Rendering, they have been used exclusively for 

35 speeding up the ray-casting rendering algorithm. The distance value are used in these systems as minimum distances 
to the next non-empty volume elements to adjust the step size along rays as they are cast into the subject. Two research 
papers that use a volumetric distance map to accelerate Volume Rendering are "Acceleration of ray-casting using 3D 
distance transforms", by K. Zuiderveld, A. Koning and M. Viergever in Proc. Visualization in Biomedical Computing, 
Chapel Hill, NC, pp. 324-335, 1 992 and "Fast surface rendering from raster data by voxel transversal using chessboard 

40 distance", by M. Sramek in Proc. IEEE Visualization '94, pp. 188-195. Washington, DC, 1994. 

[0014] Volumetric object representations have some advantages over surface-based graphical representations 
because they can incorporate complex interior structure. This is important both for visualizing internal structure and for 
modeling the physical behavior of 3D objects. However, in addition to interior structure, an accurate representation of 
surface is needed for high quality rendering, for haptic modeling, and for physics-based graphics. The subject invention 

45 addresses this need for an accurate representation of surfaces within volumetric data. 

[001 5] One reason that it is difficult to reconstruct surfaces from volumetric intensity or density values has to do with 
the discrete nature of the data. Information theory tells us that in order to reconstruct a signal exactly from sampled 
data, the signal must be sampled at twice the highest frequency component in the signal. Since a surface represents a 
sudden change in the image intensity or density value, the image signal at a surface has infinite frequency components 

so and hence no sampling rate would be high enough for exact object reconstruction. Reconstruction errors are intensified 
with gradient estimation and these errors can lead to poor surface normal estimates and significant rendering artifacts. 
[0016] The Gaussian filtering method of Wang and Kaufman prefilters the image data to reduce the high frequency 
components in a data volume and hence reduce aliasing artifacts in the gradient estimation. While this approach pro- 
duces more pleasing images than those subject to aliasing, it also results in overly smoothed and rounded surfaces. 

55 [0017] The subject invention encodes, into the volumetric data, knowledge of the presence, position and structure of 
known surfaces. Instead of storing object intensity or density values whose values usually change suddenly and dras- 
tically across a surface, the subject invention stores distances to the nearest surface element in each volumetric ele- 
ment. In one embodiment, distances are positive inside the object and negative outside of the object. 
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[0018] Note that: 

1) the gradient vector of this distance map points in the direction of the surface normal; 

2) a zero distance map value indicates the presence of a surface; and 

5 3) a positive distance value indicates that the sample is inside the object. The distance to the nearest surface point 
varies linearly in the direction norma! to a surface, and varies slowly in the direction parallel to the surface as long 
as the surface is relatively smooth, Hence, lower sampling rates and less complex gradient reconstruction filters 
can be used for accurate normal estimation with a distance map representation than with a sampled density repre- 
sentation. 

10 

[001 9] It will be appreciated that in addition to using the distance map representation for calculating surface normals 
in Volume Rendering, distance maps can also be used to calculate surface normals for haptic exploration with a force 
feedback device and for calculating object penetration distances in interaction force calculations for physically-based 
modeling. 

is [0020] As there are many sources of data for a volumetric object representation, there are many ways to calculate the 
volumetric distance map. If the source of the data is an implicit or parametric representation of the surface, such as the 
equation of a sphere or a spline patch, then either geometric reasoning, or a constrained optimization technique can be 
used to find the minimum distance from any point in the volume to the defined surface. For a polygonal surface model, 
various algorithms used in graphical collision detection can be used to calculate the distance from any point in volume 

20 to the nearest surface polygon. For a discretely sampled continuous field, such as a density field for which no further 
knowledge of the surface is known, one method to generate distance measures would be to construct a polygonal 
spline patch surface model at iso-surface values using a technique such as Marching Cubes, described by W. Lorensen 
and H. Cline, in "Marching Cubes: a high resolution 3D construction algorithm", in Proc. SIGGRAPH '87, pp. 163-169, 
1 987, and then calculate distance from this constructed surface. For a binary segmented object, the binary object could 

25 be smoothed to form a continuous sampled density field and then treated as above, or distance could be calculated 
from the binary segmented data using measurements such as the city-block distance or Euclidean distance approxima- 
tions that are described, for example, in "The image processing handbook", by J. Russ, published by CRC Press, in 
Boca Raton, FL, 1992. The resultant distance maps could also be ffltered to obtain a smoother distance field. 
[0021 ] The subject system finds application in computer graphics for generating lighting effects for volume rendering. 

30 for calculating reflected forces for volumetric objects in haptic rendering systems and for calculating physically-based 
response to interactions between volumetric graphical objects. 

[0022] With respect to generating light effects in volume rendering, the calculation of the magnitude and direction of 
light reflected from an object's surface requires an estimate of the object surface normal. The subject system for calcu- 
lating surface normals provides a means for generation these light effects for discrete volumetric objects. 

35 [0023] With respect to calculating reflected forces in haptic rendering systems, the direction and magnitude of the 
force that is reflected in haptic rendering through use of a force feedback device depends on the normal of the graphical 
object's surface at the position of the end effector of the force feedback device and the distance of the end effector from 
the object surface. The reflected forces allow one to haptically explore, or "fee", the surface of the graphical object The 
subject system of encoding distance to surface in a volumetric distance map and for calculating surface normals pro- 

40 vides a means for generating the required reflected force direction and magnitude. 

[0024] With respect to calculating physically-based responses to interactions between graphical objects, an interac- 
tion is detected when two objects penetrate each other in a physically-based simulation, and the forces caused by the 
interaction can be calculated from the depth of penetration of the two objects and the object surface normals at the con- 
tact point. The subject system of encoding distances to surfaces in a volumetric distance map and for calculating sur- 

45 face normals provides a means for detecting object penetration and for generating the required interaction force 
direction and magnitude. Object penetration of one object into a second object is detected when a point of the first 
object maps to a position in the distance map of the sec object that has a distance value which indicates that the point 
is inside the second object. The depth of penetration of the first object point into the second object can be calculated 
from the distance map of the second object. The depth of penetration can be combined with object surface normals cal- 

so culated from the distance maps to calculate the direction and magnitude of interaction forces. 

BRIEF PESCRIPTIQN QF THE DRAWINGS 

[0025] These and other features of the Subject Invention will be better understood in relation to the Detailed Descrip- 
55 tion in conjunction with the Drawings, of which: 

Figure 1 is a diagrammatic representation of a computer system used for high quality rendering; 



4 



EP0915 434 A2 



Figure 2 is a diagrammatic representation of vector components in the Phong illumination model; 

Figure 3A is a diagrammatic representation of the vector components used for calculating the lighting effects in a 
polygonal graphical object model; 

Figure 3B is a diagrammatic representation of a means for estimating the surface normal at a point, P, on a polygo- 
nal object model; 

Figure 4 is a flow chart of Volume Rendering using ray-casting with an illumination model; 
Figure 5A is a diagrammatic representation of Volume Rendering using ray-casting; 

Figure 5B is a diagrammatic representation of the vector components needed at a ray sample point, P, for calculating 
the illumination effects; 

Figure 6A is a typical data structure written in the C programming language for representing a volumetric element 
for Volume Rendering; 

Figure 6B is a data structure written in the C programming language that encodes both the element density and the 
20 distance map; 

Figure 7 is a diagrammatic representation illustrating the differences between the prior art and the subject inven- 
tion for shading in Volume Rendering, with (a) representing the two dimensional graphical object, (b) rep- 
resenting the object's distance map, (c) representing the sampled version of (a), (d) representing the 
25 sampled version of (b), (e) depicting the gradient of (c), calculated using the central difference algorithm, 

and (0 depicting the gradient of (d), calculated using the central difference algorithm; and 

Figure 8 is a block diagram of a system for generating distance maps for graphical objects from a variety of data 
sources. 
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DETAILED DESCRIPTION 



[0026] Referring now to Figure 1 , high quality rendering of a graphical object 2, is performed on a computer 4, with 
user input from the keyboard 6, and an input device 8, and displayed on a computer monitor 10. 

35 [0027] Referring to Figure 2, in computer graphics, high quality image rendering usually incorporates an illumination 
model such as the Phong illumination model which estimates the color to the light reflected at a point P from a graphical 
object as a function of the color of the object, the color of the light source, the ambient color in the room, the reflection 
coefficients of the object, the direction, L, of the light source, the direction, V, as shown, of the eye 22. and the surface 
normal, N, at P. For a given light source and a fixed viewing direction, only the object color and the surface normal vary 

40 over the surface of the object. Accurate estimation of surface normals is essential for realistic rendering. 

[0028] More particularly, it will be appreciated that the intensity of the light at each pixel of this image plane includes 
a component due to reflectance from the object surface. The refeledance intensity at point P1 is affected by the viewing 
direction V 1( from the image plane 20 to the object 22; the lighting direction, L 1( from the light source 24 to the point 
where the viewing ray intersects the object; and the surface normal, N v at that point 

45 [0029] In one embodiment, with the object being a cup, that which is necessary to accurately render the cup is to be 
able to obtain the intensity, light direction and color of a point on the cup. As to image plane 20, if one can establish a 
normal to point then one can describe the magnitude or intensity of light from that point as it relates to the viewing 
direction and the lighting direction. Since the voxel at point P<\ contains information as to color and intensity, all that is 
necessary to establish is the normal at point Pt to properly establish the intensity and color of a corresponding point on 

so the image plane. 

[0030] Referring to Figure 3A, in polygon modeling, graphical objects are represented by surface polygons, e.g. 24, 
26 and 28. When rendering polygonal models, surface illumination at a point P, 30. is calculated from the viewing direc- 
tion, V, 32, the lighting direction. L, 34, and an estimate, N, 36, of the surface normal at the point P, 30, using an illumi- 
nation model such as the Phong illumination model. Referring to Figure 3B, typically, accurate surface normals for the 
55 object are stored in the object model at the polygon vertices, 40, 42, 44, 46 and 48. Typically, the surface normal, N, 50. 
at the sample point P, 38, is estimated from the stored surface normals at the vertices of polygon 52 on which the sam- 
ple point lies. 

[0031] Referring to Figure 4, in ray casting for Volume Rendering, rays are cast into the volumetric data from each 



5 



EP0915434A2 



pixel in the view plane 54. Each ray is sampled at equally spaced intervals as it traverses the volume 56. At each sample 
point, the object color 58, transparency 60. and surface normal 62, are estimated. These are used along with an illumi- 
nation model to calculate the light reflected at each sample point 64. Once the color of each sample point is known, it 
is composited onto the ray color 66. 

5 [0032] Referring to Figure 5A. rays from the view plane, e.g. 68, 70 and 72, are cast into the volume in the viewing 
direction, sampling the color of the volumetric data at equally spaced intervals, e.g. 74 and 76 and compositing the sam- 
ple color onto the ray color. Referring to Figure 5B, in high quality Volume Rendering which uses an illumination model, 
the color of each element includes a component due to light reflected by the sample element. The color of the reflected 
light is calculated using an illumination model such as the Phong illumination model, and the resultant light color at the 

w sample point, P. 78, depends on the direction, L, 80. of the light source, the direction, V, 82, of the viewing direction, and 
the surface normal, N, 84 at P, 78. Unlike polygonal models, accurate surface normals are not usually stored in the vol- 
umetric data. Hence, the surface normal, N, 78. at the sample point must be calculated from the density values stored 
in data elements near the sample point. Referring to Figure 6A, typically, in Volume Rendering, each element in the vol- 
ume has a structure such as the one shown here, written in the C programming language. Each element contains a 

is value 86 for the object density or intensity. Referring to Figure 6B. in the subject invention, the data structure for each 
element has an additional field 88 so that a value for the distance to the nearest surface can be stored for each element. 
In the subject invention, this distance is used to detect the presence of a surface and to estimate its normal direction. 
In one embodiment, the distance is positive inside of the object negative outside of the object and zero on the object 
surface. The gradient of this distance field points in the direction of the surface normal. Surface are detected where the 

20 magnitude of the distance is below a small threshold. 

[0033] Referring to Figure 7, differences between the prior art for estimating surface normals in Volume Rendering 
and the method of the subject invention are illustrated in 2D. 7(a) is a 2D object. The object is present where the image 
is black as illustrated at 90, and absent where it is white, as illustrated at 92. Figure 7(b) is the continuous distance map 
that corresponds to the object in Figure 7(a). Note that it has a value of 0.5 on the surface, or edge, of the 2D object 94; 

25 values that increase with distance from the surface and ranging from 0.5 to 1 inside the object; and values that decrease 
with distance from the surface and ranging from 0.5 to 0 outside of the object. Figure 7(c) is the sampled version of Fig- 
ure 7(a), where sample points are black if they fall inside the object, e.g. 96, and white if the fail outside of the object, 
e.g. 98. Figure 7(d) is the sampled version of Figure 7(b), where the value of each sample point is the color of the dis- 
tance map at the same position. Notice that the intensity values of the sample points in Figure 7(d) vary smoothly 

30 across the image. Because of this smooth variation, values of the distance field that lie between sample points can be 
reconstructed more accurately and with a lower order interpolation filter than values of the sampled image density field 
of Figure 7(c) can be reconstructed. Figure 7(e) is a sampled estimate of the gradient field of Figure 7(a) that was cal- 
culated from the sampled image of Figure 7(c) using a central difference gradient estimator. Note that the direction of 
the sample gradients can vary significantly for neighboring samples, e.g. 98 and 100. or 102 and 104. Using this gradi- 

35 ent estimate in Volume Rendering would cause significant artifacts in the image. Figure 7(f) is an estimate of the surface 
normal directions of Figure 7(a) calculated using a central difference gradient estimator from the sampled image of Fig- 
ure 7(d). The direction of the surface estimate is consistent throughout Figure 7(f) and it accurately represents the nor- 
mal of the 2D object surface. TTie magnitude of the surface normal can be derived from the distance map in Figure 7(d). 
In this example, the surface normal magnitude is large when the distance map intensity is close to 0.5 and small other- 

40 wise. 

[0034] Referring to Figure 8. there are different ways to calculate an object's distance map 106. depending on the 
object representation 108. If the object is represented implicitly or parametricaJly as at 1 10, such as the equation of a 
sphere or a spline patch, then either geometric reasoning or constrained optimization 1 12 can be used to find the min- 
imum distance from any point in the volume to the defined surface. If the object is represented by a polygonal surface 

45 model 1 1 4, various algorithms used in graphical collision detection can be used to calculate the distance from any point 
in volume to the nearest surface polygon 1 16. For a discretely sampled continuous density field 118, one method to 
generate distance measures would be to construct a polygonal surface model at iso-surface values 120, and then cal- 
culate distances from this constructed surface 116. For a binary segmented object, the binary object could be 
smoothed as illustrated at 124 to form a continuous sampled density field and then treated as above as illustrated at 

so 120 and 1 16, or Euclidean distances could be calculated from binary segmented data 126. The resultant distance maps 
could also be filtered to obtain a smoother distance field 1 28. 

[0035] A program listing in C is presented hereinafter to illustrate the utilization of the subject system for use in gen- 
erating lighting effects for volume rendering: 

55 
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/* volRender.c 

* 

* Basic volume rendering program. Takes a data volume consisting of an 

* intensity and a distence-to*surf ace measure for each voxel and projects 

* the volume onto a 2D image. 

* Revised from earlier projection programs. 
« 

* Calculates object color, and transparency using the voxel intensity values 

* and shading based on a central difference of the distance values. 

* Compositing is done in front-to-back order and terminates at a max opacity. 

* July. 1997 by S. Gibson. 

* MERL Cambridge. 
* 

•/ 

/• includes •/ 
flinclude "volRender .h" 

/• Vis_Display global variables */ 

static int xSize. ySize. zSize; 

static int imgSize; 

static int nLights; 

static int vertexOf f set ( 8 1 ; . 

static float voxelSizeX. voxelSizeY. voxelSizeZ; 
static float scaleJC. sceleY, scaleZ; 
static float R{3}{31, T(3); 

static float kd. ks, specularExp. visibility; 
static float imageGain; 

static VisUght lights (MAXLIGHTS] ; 
static VisFVector view; 



30 * main 



35 



40 



void main {int argc, char *"argv) 



( 



char da ta Fi leName (60); 
char clutFi leName ( 60 ] ; 
char tlutFileName ( 60 ] ; 



int i, j. 
int size, 
int fd, n; 



k , index ; 
imgArea ; 



XEvent event; 



45 



50 



VisVoxel * voxel; 
VicColor *ch«rlmnge; 
VisFColbr -imager 
VisColorLut clut; 
VisTranspLut tlut; 



/* 3D data volume •/ 

/' projected image (char) •/ 

/• projected image (floating point) ■/ 

/• array of. color lookup tables •/ 

/* array of transparency lookup tables •/ 



/* initialize the rendering parameters V 

InitProjection(argvU) . data? ileName. clutrileName . tlut Filename) ; 

/* allocate space for the cnta and image arrays •/ 
size s xisize * ySize * zSize; 

voxel * (VisVoxel *! malloc(size * sizeof (VisVoxel) ); ' 



55 
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imgArea « ingSize * imgSize; 

image s (VisFColor *) malloc (imgArea ■• sizeof (VisFColor ) ) ; 
charlmage - .(VisColor *) malloc (imgArea * sizeof (VisColor) ) ; 

/* opeu and read the data file •/ 
fd = open(dataFileName, OJIDONLY) ; 
read(fd. voxel, size * sizeof (VisVoxel) ) ; 
close(fd) ; - . 

/« clear the imafle to block # / 
VisClear Image (image ) ; 

/* render the volume data •/ 
VisProject (voxel . clut, tlut. image); 

/* scale the colors of the image to chars, correct for overflow •/ 
ColorsToChar (image . char Image , imgArea ) ; 

/• write the image to a file */ 
fd * creat ("image. raw' . PhOD£) ; 

writelfd, charlmage. imgSize * imgSize # si2eof (VisColor) ) ; 
close (fd) ; 

/* free the data allocation •/ 
free( (VisVoxel •) voxel); 
free ( (VisFColor *) image); 
free ( (VisColor * ) char Image) ; 



VisProject — 

This procedure calculates the projected image by casing rays from the 
image plane through t the object, calculating the gradient, color, 
transparency, and shaded color contribution from each sample point 
along the ray and then compositing these values in a front- to- bach 
ordering. 



void , 
VisProject (VisVoxel * voxel. VicColorLut clut. Vis-TranspLut tlut. 
VisFColor • image) 

( 

int i. j. k. ii. jj, kk; 
int index, nlndex; 

int planeSize; sj 
int neighbor Index (8) , values (8); jj 

float x. y, z; 

float dx, dy. dz. ddx. ddy. ddz; 
float tempo, tempi, temp2. temp3 ; 
float w(B) ; 

float transp. opacity. accTransp; 
float mag. sign, dist; 

VisFColor color. shadingColor ; 
VisFVector grad; 

planeSize = xSize * ySize; 

/• 1) calculate contributions to the image plane by sampling along a ray 
• through the data in the kk direction in front- to-back order 
-/ 
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2) for each pixel (ii.jj) in the image plane, 1) find the 

corresponding coordinates (x.y.z) in the data volume for the kkth 
sample point along the ray, 2) trilinearly interpolate attributes 
of the sample point from the surrounding neighbors. 3) accumulate 
the color along the ray by compositing the color of the sample 
point onto the image plane pixel. 

ndex = 0; 

or (jj * 0; jj < imgSize; jj*+) ( 
for (ii = 0; ii < iwgSize; ii*+) .( 
accTransp = 1.0; 
kk = 0; 

color. r = color. g * color. b = 0.0; 
transp = 1.0; 

while ((kk < imgSize) (accTransp > 0.001)) ( 

/• 3) calculate the (x.y,2) coordinates of (ii.jj,kk) in image space •/ 
x = Rt0|(0j * ii + R(0J(1] 4 jj + n[0](2] • kk ♦ T(0J; 
y = R { 1 1 1 0 ) * ii * R[1][1J * jj + 1UDI2J • kk ♦ T(1J; 
2 = R[2){0] * ii * RC2U1) * jj + n|2J(2J • kk + T(2] ; 

/• make sure the (x.y.z) lies inside of the data volume */ 
if ((x >= 0) (x < xSize-1) Uk 

(y >= 0) t« (y < ySize-l) £t 

(z >= 0) Li. (z < zSize-1) ) [ 

/* 4) calculate coordinate dependent weights and constants •/ 

i = (int) x; j = (int) y; k = (int) z; 

dx = x-i ; dy = y- j ; dz = z-k; 

ddx = 1-dx; ddy = 1-dy; ddz = 1-dz; 

tempo = ddx'ddy; tempi = ddx'dy; temp2 « dx*ddy; temp3 = dx<dy; 

/* weights for trilinear interpolation where the indices correspond 

♦ to the neighbors in the following way: D->000. l-xD01, 2->010. 

• 3->031, 4->100. 5->101, 6->110. and 7->lll. 
•/ 

w(0) = tempO-ddz; w(l) = tempO'dz; w(2J = templ*ddz; w[3J = templ'dz; 
w{4) = temp2*ddz; w(5) = temp2«dz; w£6J = temp3'ddz; w(7) = temp3*dz; 

r 

nei ghbor Index ( 0) = i + j'xSize + k*planeSi2e; 

dist - VisCetSampleDist (voxel , w, neighborlndex ( 0 J ) ; 

if (dist >= -0) ( /* odd color only for interior elements •/ 

/* indices of B surrounding voxels where the indices correspond to 

• the neighbors in the following wayu 0->000. 1->001, 2->010, . 

• 3->011. 4->100. 5->101. 6->110, and 7->Ul. 
•/ * 

for (nlndex = 1; nlndex < 8; nlndex »•+) ( 
neighborlndex I nlndexj = 

neighborlndex! 0] + vertexOf I set tnlndexl ; 

) 

for (nlndex = 0; nlndex < B; nlndex«-+) ( 

values [nlndex) = (int) voxel (neighborlndex (nlndex) ) .value; 

J 

/* 

• 5) calculate the gradient, color, transparency, and shaded 

• color of the sample point from the surrounding neighbors. 
V 

VisGetSampleColor (values, w, ciut. tcolor) ; 
transp = VisGetSarapleTransp (values, w. tlut); 
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if (nLights > D) ( 

it 0»ag = VisGetSampleCrad (voxel, w; i, j, k. igrad) >= 0} ( 
. /• find direction of grod relative to the view direction ^ 

• If positive, this sample point contributes to the shading 

*/ 

sign = VectorDot (grad, view); 
if (sign >= 0) ( 

shadingColor = color; 

VisGetShadedColor (grad. mag. ^shadingColor ) ; 

color. r += shadingColor .r; 
color. g += shadingColor .g; 
color. b shadingColor .b; 

) 

) 

) 

/* 6) project the voxel onto the view plane */ 
impge[indexl .r += accTransp • color.r; 
image ( index) .g ♦ = accTransp * color. g; 
imoge ( index) .b accTransp • color. b; 
accTransp *= transp; 

) 

) 

; 

) 

index**; 

) 

) 

return; 

) 



VisGetSajnpleColor 

This procedure gets the sample color by trilinearly interpolating the 
colors of its 8 nearest neighbors after they. have been determined 
using a color lookup table. 



V 
void 

VisGetSnrapleColor (int *values, float 'weight, VisColorLut clut. 
VisrColor *color) S 
( * 
int i; t 
VisColor temp; 

color->r = color - > g = color->b = 0; 

for (i = 0; i < 8; i++)' ( 

temp = clut. color (values (i] ) ; 

color->r 4= visibility * weiphtUJ * (float) temp.r / 255.0; 
color->g visibility • weight [i) * (float) terop.g / 255.0; 
color->b += visibility • weight (i 1 " (float) temp.b / 255.0; 

J 

return; 

) 
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r 



* VisGetSaapleTransp 

* This procedure gets the cample transparency : by trilineaxly interpolating 

* the transparencies of its B nearest neighbors after they hove been 
determined using a lookup table. 



10 

float 

VisGetSarapleTransptint 'values, float 'weight. VisTranspLut tlut) 
( 

int i ; 

float Lrensp = 0; 

for U » 0; i < B; i++) { 

transp ' weight (i I * tlut . transp (values [ i J I ; 

) 

transp /* 255.0; /* »»x char value is 255 */ 

return transp; 

20 J 



15 



• VisGetSampleDist -- 
25 * 

* This procedure gets the sample distance by trilinearly interpolating 
the distances to near-st surfaces of its B 'nearest neighbors. 



30 float 

VisGetSampleDist {Vis Voxel 'voxel, float -weight, int index) 

C . 
int i; 

float dist = 0; 

35 for {i » 0; i < li; i++) 

dist += weight (i) ' voxel (index ♦ vertexOff set [i} } .dist; 

■ i 

return dist; 

} 



• VisGetSampleGrad — 



45 



50 



This procedure gets the sample gradient from the distances of the 
* neighboring elements using a central difference method. 



float 

VisGetSampleGradtVisVoxel 'voxel, float 'weight, 

int voxell, int voxelJ. int voxclK. VisFVector *grad) 

I 

int i, j, k. index; 
float mag; 



55 
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/• apply central difference operator to trilinearly interpolated distances 

of the voxel's 6 nearest neighbors to estmate the image gradient */ 
/* note: samples must lie inside the data volume */ 
if ((voxell > 1) kk (voxell < xSi2e-2) kk 
(voxelJ > 1) kk (voxelJ < ySize-2) kk 
(voxelK > 1) kk (voxelK < 2Size-2)| ( 
index = voxell ♦ voxelJ*xSize ♦ voxelK*xSize*ySi*e; 

grnd T >x s - VisGetSampleDist (voxel . weight, index ♦ 1) 4 

v'ifiCetSompleUist {voxel, weight, index - 1); 
grad->y * - VisGetSampleDist (voxel . weight, index ♦ xSize) ♦ 

VisGetSainpieDist (voxel , weight, index - xSize); 
grad->2 = - VisGetSampleDist {voxel , weight, index «• xSi2e"ySize) + 

VisGetSampleDist (voxel, weight, index - xSize*ySize) ; 

mag = sqrt(grad->x * grod->x * grad->y * grad->y + grad->x * grad->z) ; 

) 

else ( 
mag s 0; 

) 

if (mag > 0) { 
gred->x /= mag; 
grad->y /= may; 
gritd->2 /= mag; 

) 

else ( 
mag = 0 ; 

grad->x = grad->y = yrad->z = 0; 

J 

return mag; 



VisGetShadedColor — 

This procedure calculates shaded color using a Phong shading algorithm 
with diffuse and specular reflection. 



void VisGetShadedColor (VisFVector grad, float magi. VisFColor • color) 

int i; ji 
float dScale, temp, sScaie; 
VisFColor shadeColor; 

shadeColor. r = shadeColor.g = shadeColor .b = 0; 

if (visibility !« 0) ( 

for (i s 0; i < nLighcs; i+«-) ( 

/* get the diffuse lighting scaling factor */ 
temp ■ VectorDot (grad. lights (i J .lightVector J ; 
if (temp > 0) 

dScale = mag * kd * temp; 
else 

dScale = 0; 

/" get the specular light scalinu Cactor. Note the use of the halfway 
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* vector assumes that the light and the viewing position are infinitely 

* far away (see Foley, Van Dam, et al.). 
+ / 

temp = VectorDot (grod. lights (i | .hwVector) ; 
•if (temp > 0) 

sScale = mag * ks • pow(temp. speculnrExpJ ; 
else 

sScele = 0; 



10 



shedeColor.r *s visibility * lights(i) 
shodeColor.g visibility * lights (i) 
shodeColor . b -»•= visibility ** lights (i J 



r «\(dScale 



(dScale 
(dScale 



color->r sScalc) 
color->g + sScale) 
color->b sScale) 
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•color s shadeColor; 
return; 



20 [0036] Having above indicated several embodiments of the Subject Invention, it will occur to those skilled in the art 
that modifications and alternatives can be practiced within the spirit of the invention. It is accordingly intended to define 
the scope of the invention only as indicated in the following claims. 

Claims 

25 

1 . In a volumetric data representation system for graphical objects in which graphical objects are generated from data 
elements, 

means for encoding the distance of a data element to the nearest object surface point; and, 

30 

means responsive to the distance of said data element to the nearest object surface point for estimating a sur- 
face normal thereof. 

2. The system of Claim 1 , and further including means responsive to the distance of said data element to the nearest 
35 object surface point for estimating the magnitude of said surface normal. 

3. The system of Claim 1 . wherein said distance encoding means includes means for generating a discrete volumetric 
distance map and means for generating a gradient vector at any point in the space that contains said surface from 
said distance map, said gradient pointing in the direction of a surface normal. 

40 

4. The system of Claim 3, wherein a zero distance map value indicates the presence of a surface. 

5. The system of Claim 3. wherein a positive distance value indicates that said data element is inside said object. 

45 



50 



55 
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Fig. 3A 
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Fig. 3B 
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FOR EACH PIXEL IN THE VIEWING PLANE 



CAST RAYS INTO THE VOLUME 
54 



SAMPLE VOLUMETRIC DATA 
AT INTERVALS ALONG THE RAY 
56 



AT EACH SAMPLE. 



50- 



INTERPOLATE OBJECT COLOR 



60~ INTERPOLATE OBJECT TRANSPARENCY 



62- 



ESTIMATE OBJECT SURFACE ' 



64- 



USE SURFACE INFORMATION 
MODEL TO ESTIMATE THE 
• REFLECTED LIGHT 



66- COMPOSITE SAMPLE RESULT ON RAY COLOR 



Fig. 4 
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TYPE DEF VOLUME ELEMENT STRUCT { 

06 - INT DENSITY; 

} VOLUME ELEMENT; 



Fig. 6A 



""TYPE DEF VOLUME ELEMENT STRUCT { 
INT DENSITY; 

88 INT DISTANCE; 

} VOLUME ELEMENT; 



Fig. 6B 



20 



EP0915434 A2 




Fig. 7A Fig. 7B 
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Fig. 7C 
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Fig. 7D 




Fig. 7E 
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Fig. 7F 
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OBJECT REPRESENTATION DISTANCE CALCULATION 



DISTANCE MAP 
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IMPLICIT/ 
PARAMETRIC 
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GEOMETRY/ 
CONSTRAINED OPTIMIZATION 
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POLYAGONAL 
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DISTANCE TO CL05EST POLYGON 
(FROM PHYSICS-BASED GRAPHlCSI 
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SAMPLING 
CONTINUES 
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BINARY 
SAMPLING 



ISO SURFACE 
GENERATION 
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SMOOTHING 
FILLER 
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MAP 
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Fig. 8 
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